%% The actual Travel model and modifications
%
% Input:  PreferenceFactor (used for minimization),cfg,exvars,depvars,indepvars,r,t
% Output: Residuals (depeding on calibration type; only during calibrations)
%         depvars structure
%
function varargout = ModelSimple(cfg,exvars,depvars,indepvars,r,time,varargin)
%% Destinguish between calibration or validation and specific type of run
for n = 2:2:length(varargin)
    eval( [varargin{n-1} '(r,1,:)=varargin{n};']);
end
modes  =  indepvars.UseModes(r,1,:); % get rid of modes not used by Schafer (e.g. High speed train)
indepvars.signES(1,1,1:4) = [-1,-1,-1,-1];
indepvars.signEY(1,1,1:4) = [-1,1,1,1];

for t=time
    if t==1
        depvars.TravelDemand(r,t,modes) = exvars.TravelDemand(r,t,modes);
        
    else
        depvars.TravelDemand(r,t,modes) = depvars.TravelDemand(r,t-1,modes) ...
            .*(exvars.Income(r,t)./exvars.Income(r,t-1)).^(indepvars.signEY(1,1,modes).*indepvars.IncomeElasticity(r,1,modes))...
            .*(exvars.TravelPrice(r,t,modes)./exvars.TravelPrice(r,t-1,:)).^(indepvars.signES(1,1,modes).*(indepvars.PriceElasticity(r,1,modes)));
    end
    
    depvars.ModeSplit(r,t,modes)         ...
        = depvars.TravelDemand(r,t,modes) /sum(depvars.TravelDemand(r,t,modes),3);
%% Errors
    if t<36 % ensures function is only called with historic data
        depvars.ModeWeight(r,t,modes) = sqrt(exvars.TravelDemand(r,t,modes))./sum(sqrt(exvars.TravelDemand(r,t,modes)),3);% Ensures not all modes are equal, but that car is not the only important mode
        depvars.ObjectiveFunction(r,t) = sqrt(sum((depvars.TravelDemand(r,t,modes) - exvars.TravelDemand(r,t,modes)).^2 ...
        ./ exvars.TravelDemand(r,t,modes).^2 .* depvars.ModeWeight(r,t,modes).^2,3)); % used for optimization
        depvars.Resid(r,t)             =  depvars.ObjectiveFunction(r,t);
    end
end

%% Write output
    varargout{1} = mean(depvars.ObjectiveFunction(r,time));
if nargout>=2
    varargout{2} = depvars;
end
if nargout >=3
    varargout{3} = indepvars;
end

%% Visualizes minimization process
if cfg.LiveMinimizeVisuals{1} > 1 && ismember(time(1),cfg.LiveMinimizeVisuals{2}); % bit clumsy to make the code faster
    GeneralPlotting('LiveMinimizeVisuals','Dont save',cfg,exvars,depvars,indepvars,r,time);
elseif cfg.PlotMinimize == 1
    GeneralPlotting('SaveToGlobal',' ',cfg,exvars,depvars,indepvars,r,time);
end
    
end